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ABSTRACT 



Many scenarios have been proposed for the origin of the supermassive black holes 
(SMBHs) that are found in the centres of most galaxies. Many of these formation 
scenarios predict a high-redshift population of intermediate-mass black holes (IMBHs), 
with masses M, in the range IO^Mq < M, < W^Mq. A powerful way to observe these 
IMBHs is via gravitational waves the black holes emit as they merge. The statistics of 
the observed black hole population should, in principle, allow us to discriminate between 
competing astrophysical scenarios for the origin and formation of SMBHs. However, 
gravitational wave detectors such as LISA will not be able to detect all such mergers nor 
assign precise black hole parameters to the merger, due to weak gravitational wave signal 
strengths. In order to use LISA observations to infer the statistics of the underlying 
population, these errors must be taken into account. We describe here a method for 
folding the LISA gravitational wave parameter error estimates into an 'error kernel' 
designed for use at the population model level. The effects of this error function are 
demonstrated by applying it to several recent models of black hole mergers, and some 
tentative conclusions are made about LISA's ability to test scenarios of the origin and 
formation of supermassive black holes. 

Subject headings: Black Hole Physics - Early Universe - Gravitational Waves - Methods: 
Statistical 



Introduction 



There is now substantial evidence (e.g., Kormendy & Richstone 1995[ [Richstone et al. 1998 



Bender|[2005 Rees||2002 20031 for the existence of supermassive black holes (SMBHs) in the nuclei 



of most galaxies, the black hole in our own galaxy being the best studied and most clearly justified 
of these objects. However, the origin of these black holes remains an unsettled question. In one 
scenario, the more massive black holes formed from the merger and coalescence of smaller 'seed' 
black holes that were created in the very early Universe (e.g., Madau Sz ReespOOl l. Several models 
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utilizing this process have been proposed and numerically simulated (e.g., 


Haehnelt &: Kauffmann 


2000 


Volonteri et al. 


2003 


I. In (e.g., 


Volonteri et al. 


2003 


Tanaka &: Haiman 


2009 


1, typical seed 



black holes are the remnants of Population III stars with masses m, ~ 100 — 3OOM0 formed at 
high redshift (e.g., z ~ 20 - 50). Thus, these models predict an evolving population of intermediate 
mass black holes (IMBHs), with masses between ~ 100 to lO^M0j^In another version (Begelman 
et al.|[2006 ) , the seed black holes are formed due to direct collapse of the cores of pregalactic halos 
through a 'quasi-star' stage, resulting in a more massive seed population (~ 10^ - IO^Mq). This 
scenario predicts far fewer IMBHs in the early universe. 

Black holes in the IMBH mass range are extremely difficult to detect with the usual electromag- 
netic observation techniques, making it very difficult to verify a particular formation and evolution 
scenario and, especially, to discriminate between various models. However, the mergers themselves 
produce gravitational waves in the Low Frequency (LF) band, from 10~^ Hz to 10~^ Hz, probed by 



the proposed LISA mission (Jennrich 2004 1 . Observations of the amplitude, frequency chirp, and 
harmonic structure of the gravitational wave waveform enable both the luminosity distance and the 
individual masses of the black holes in the binary system to be determined. LF gravitational wave 
observations thus provide a probe of the cosmological spectrum of black holes, and allow tests of the 
population models to be made. In this paper, we investigate how gravitational wave observations 
of coalescing massive black hole binaries may be used to discriminate between models of massive 
black hole populations and determine the merger history that has led to the observed population 
of SMBHs. 

We thus take a slightly different direction, compared to other recent works investigating LISA 
detections of black hole coalescences. Most of these take the somewhat speculative black hole popu- 
lation models and calculate the number of coalescences that each model predicts would be observed 
by LISA. The scientific goal of these papers has clearly been the observation of the coalescence 
itself. In this paper, we turn the problem around and ask, 'What might LISA observations of 
binary MBH coalescences tell us about the otherwise uncertain population models?' 

The organization of the paper is as follows. In Section 2, we discuss some of the population 
models and the motivation behind them. Section 3 describes the parameters relevant to detection of 
a black hole binary and draws a distinction between 'population' parameters, which are relevant to 
the population models, and astrophysically uninteresting 'sample' parameters, which vary randomly 
for each sample drawn from a population. We also review the response of the LISA gravitational 
wave detector to massive black hole binary coalescences, with emphasis on the ability of the detector 
to determine the parameters of the binary from the gravitational wave waveform. The bulk of our 
work is presented in Section 4, which describes a Monte Carlo calculation of an 'error kernel', which 
is marginalised over the sample parameters. This error kernel, K{\i, Aj), is the average conditional 
probability that a source will be detected with population parameters \i in the LISA data, given 



Since our studies include both IMBHs, with mass ~ lOOAf0 to 10^ Mq, and SMBHs, with mass ~ lO''Af0 to 
10® Mq, we adopt the terminology 'massive black hole' (MBH) to cover the entire rang e from ~ lOOA/© to IQPMq. 
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the existence of an astrophysical source with population parameters Aj. The marginahsation over 
sample parameters sets this apart from previous work, and the resulting error kernels can be applied 
directly to model coalescence rates, producing a new set of coalescence rates as functions of the 
best-fitting, 'detected', parameter values. Finally, in Section 5, we take several population models 
from the literature and discuss how the error kernel may be used to produce a measure of how well 
the models may be discriminated from one another. Although the models we use are incomplete, 
being based only on summaries available in the literature, we nevertheless draw a few tentative 
conclusions about LISA's ultimate ability to distinguish between black hole population models in 
Section 6. 



2. Astrophysical Populations of Black Hole Binary Systems 



Observations of the high redshift quasar population (Fan et al. 2001 Stern et al.||2000 Zheng 
et al.|[2000 Becker et al.||2001 1 suggest that a population of SMBHs has existed since early epochs 
~ 6). The local census of SMBHs has been increasing in recent years (Tremaine et al. 20021, 



driven by a growing body of observational evidence linking the mass of SMBHs with observational 
properties of their host galaxies. Early studies revealed a rough correlation between SMBH mass 
and the bulge luminosity of the host galaxy ( Kormendy &: Richstone||l995 Magorrian eraI][T998| ). 
A much stronger correlation was later discovered between the SMBH mass and the stellar velocity 



dispersion in the galactic core, the so-called 'M-cr' relation (Gebhardt et al. 2000; Ferrarese Sz 



MerrittpOOO 


Tremaine et al.||2002|. 


2001 Tremaine et al. 


2002 


1 gives the 



(1) 



The observational data supporting the M — a relation currently spans a mass range from ~ lO^M© 
to ~ IO^Mq. 

Given this observational evidence for the existence of a SMBH population, the question arises: 



how did these objects come to be? Several scenarios are proposed (see Djorgovski et al. 20081: 



1. direct gravitational core collapse of pregalactic dark halos, 

2. growth from seed black holes through merging and coalescences over time, 

3. gravitational runaway collapse of dense star clusters, 

4. primordial BH remnants from the big bang. 

In case [T] SMBHs can form very early in the Universe through direct collapse of dark matter 
halo with mass of ~ IO^Mq or larger (Bromm Sz Loeb 2002 1. In case|2] on the other hand, (Madau Sz 



Reesl[200T| [Haehnelt Kauttmaniil[2000l |Volonteri et~aL][2003| [Tanaka Haiman|[2009l |Begelman 



et al. 20061, seed black holes produced in the early universe are significantly smaller and grow 



by accretion, coalescences and merging, leading to the population of SMBHs seen in the Universe 
today. In most cases (e.g., Volonteri et al.||2003 | the seed black holes have mass less than ~ SOOM©, 
and are the remnants of ordinary Population III stars. Detailed stellar evolution calculations have 
recently found, however, that these seed black holes can be remnants of very massive Population 



III stars, referred to as CVMSs (Ohkubo et al. 2006 2008: Tsuruta et al. 2007 Ohkubo et al. 



2009 Umeda et al. 20091. These metal-free CVMSs evolve quickly and then collapse in the early 
universe, yielding an IMBH population with masses in the range of ~ SOOM© — 1O,OOOM0. In 



another version (Begelman et al. 20061, even more massive black holes can form from direct core 



collapse of mini halos through a quasi-star phase, instead of through ordinary stellar evolution, 
leading to a population with masses in the range ~ IO^Mq to ~ IO^Mq. 

The scenario involving direct collapse is not difficult to distinguish from that involving Popu- 
lation III stars, because the former predicts only a small number of the more massive IMBHs while 
the latter predicts considerably more IMBHs with a wider range of masses and redshifts. In case 
ll (e.g., [Ebisuzaki et~aL]|2001| [Portegies Zwart et al.|[2004a||bl ), IMBHs can be formed at any time 
in dense star clusters and grow by merging in the given environment. Such a process can produce 
a low level population of mergers at all redshifts. Case |4] is theoretically possible, but currently it 
will be hard to test by observations. 

In this paper, we concentrate on versions of case [2] although the techniques used are general 
and should also apply to discrimination between other MBH population scenarios. Specifically, we 



use the results published in Sesana et al. (20071, which give coalescence rates for models discussed 
in I Volonteri et al.l (|2003|) , iKoushiappas et al.l (|2004|) , and iBegelman et al.l (|2006|) . In these models. 



the evolution of the seed population through time is carried out numerically ( Somerville & Kolatt 



1999[ [Volonteri et al.||2003[ ). In [Volonteri et al.| ( [2003l ), for instance, the evolution of a SMBH's host 
galaxy halo is first computed, starting at the present day and working backwards in time to z ~ 20. 
Then, the progenitor haloes are seeded with BHs of ISOMq, and the BH merger history is traced 
forward in time through the halo merger history. Major differences between different versions of this 
scenario come from the varying astrophysical assumptions employed - e.g., the initial conditions, 
environment, age (the redshifts when the seeds were produced), seed mass, dynamics, hierarchial 
growth, etc., built into the process. 



3. Gravitational vi^ave binaries 

3.1. Binary Parameters 

Each black hole binary may be characterized by a number of parameters, which are usually 
divided into two categories. The first category is the intrinsic parameters which have to do with 
the local properties of the binary in its rest frame. They are mi, the mass of the primary, m2, 
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the mass of the secondary, and either the initial orbital separation, a, frequency, / (related to a 
by Kepler's third law) or time to coalescence, ic; (related to a by the quadrupole formula)]^ A set 
of mass parameters equivalent to mi and m2 but more directly related to the gravitational wave 
waveform, is the chirp mass Mc = {jni'm2)^/^{jni + m2)~'^^^ and the symmetric reduced mass ratio 
rj = (17111712) /{mi + 7712)^. Note that, in the gravitational wave signal, the masses are scaled by 
the redshift, so that the natural mass variables for gravitational wave data analysis are redshifted. 
For instance, the redshifted chirp mass is Aic x (1 + -2), although the reduced mass ratio remains 
unchanged since it is dimensionless. The remaining parameters fall into the second category, and 
are called extrinsic parameters. These have to do with the binary's location and orientation with 
respect to the LISA constellation. They are the luminosity distance Di (or equivalently, the redshift 
z), the principle gravitational wave polarization angle ip, the binary inclination t, the sky location 
angles 9 and 0, and the initial phase of the binary orbit ^q. 

Redshift and luminosity distance are used interchangeably as source parameters, with the 
relationship between them determined from the standard WMAP cosmology (Om = 0.27, ilvac = 
0.73, Orad = 0.0, and a Hubble constant of 71 km/s/Mpc). 

Since the predictions of the population studies are given as functions of masses and redshift, 
the binary parameters are best divided into two sets in a different way, for purposes of this paper. 
The first set, consisting of ^Ac■, i], and z, are what we will call population parameters, since these 
are the parameters that characterise the population model predictions. The remaining parameters, 
ip, i, 6, (j), tc, and ^0, represent particular samples drawn from the population model and will 
be referred to as sample parameters. Sample parameters have distributions that are essentially 
stochastic and contain no useful information about the astrophysical processes which give rise to 
the black hole population. 

Despite the fact that the sample parameters are not part of the intrinsic astrophysical model, 
they can have a dramatic effect on LISA's source characterisation capabilities because all parameters 
must be fit to the data in the process of extracting the population parameters of interest. A 
key component of this analysis is therefore to calculate average LISA error distributions for the 
population parameters by averaging over a Monte Carlo ensemble of many sources, each having 
randomly-chosen values for the sample parameters. Such error distributions are referred to as 'error 
kernels'. 



3.2. Gravitational Waves from Binary Systems 



Calculation of the detectability of binary systems via gravitational wave emission is a standard 
problem in the gravitational wave community; see, for instance, Flanagan Sz Hughes ( 1998a|b ), and 



■^In general, the list of intrinsic parameters also includes the black hole spins and orbital eccentricity. Here, 
however, attention is restricted to the simplified case of non-spinning black holes in circular orbits. 
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Cutler &; Flanagan (1994 1 . We review the problem here for convenience and locality of reference. 



Spin interactions and higher order post-Newtonian corrections are neglected in this paper, 



although they can have an appreciable effect on parameter error estimation. For instance, Hellings 



Sz Moore ( 2003 1 have shown that the inclusion of higher harmonics of the waveform would improve 



the determination of the mass parameter rj, while Lang Sz Hughes (2006) have shown that spinning 



black holes produce a modulation in the signal that leads to a tighter bound on the sky position of 
the binary. 

The dimensionless gravitational wave strain produced by a circularised binary can be written 
as the superposition of two independent polarization states /i+ and hx ■ The polarization amplitudes 
can be expanded in terms of harmonics as 



(2) 



where t = t — k ■ x locates the surface of constant phase for a gravitational wave propagating in a 
direction k, and $(t) is the phase of the binary orbit, as observed at the LISA detector. When the 
binary is far from coalescence, the dominant emission is the n = 2 quadrupole, which is 



hxir) 



2/3 



cos i sin [2$(r)] 



Dl 
Dl 



(3) 



where / is the fundamental quadrupole frequency, / = 2[d(^ / dt) / {2tt) . We note that /i+ and 
are still functions of the observed time r. 

The response of the LISA detector to the two polarizations of a gravitational wave from a 
binary is given by 



y(T) = F+{e, 0, i, iP, r)/i+(r) + Fx (6, cj), l, V, t)/ix (r) 



(4) 



where and Fx are the LISA form factors that depend on the position and orientation of the 
source relative to the time-dependent LISA configuration. 



3.3. Detecting Black Hole Binaries 

One measure of the ability of the LISA detector to observe a binary signal is the signal-to-noise 
ratio, defined as 



(SNR)^ 



\h{f)\' 
5'lisa(/) 



df 



(5) 
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where = + |/ix(/)P, with h^{f) and hx{f) being the Fourier transforms of the 

polarization amphtudes in equation |3j and where >S'lisa(/) is the apparent noise level of LISA's 



Standard Curve Generator (Larson (2000), hereafter SCG), an estimate that averages the LISA 
response over the entire sky and over all polarization states and divides the LISA instrument noise, 
Snif), by this averaged response. 



Previous treatment of LISA observations of binary black hole populations ( Sesana et al.|[2007 | 



have employed this measure of detectability, while others (Sesana et al. 2004 1 have used a charac- 



teristic strain he, following the prescription of Thorne (1987 1 . In this measure, the raw strain /i of a 



source is multiplied by the average number of cycles of radiation emitted over a frequency interval 
A/ = l/Toi)s centred at frequency /. The amplitude of the characteristic strain is then compared 
directly against the 1-year averaged strain sensitivity curve from the SCG to produce an SNR. In 
either case, a source is considered detectable if the resulting SNR exceeds some standard threshold 
value (typically between 5 and 10). 

While interesting for planning LISA data analysis pipelines, these SNR estimates fail to address 
the fact that a detection is of little use for comparison with astrophysical theory if the parameters 
of the binary are poorly determined. In particular, unless the masses and redshifts of the detected 
black holes are measured, the observations cannot be compared with the black-hole evolution mod- 
els. A more complete analysis that incorporates the effects of uncertainty in the binary parameters 
is required. 

Parameter error estimation for black hole binaries detected via gravitational wave emission has 



been discussed by many researchers (Cutler Sz Flanagan 1994 Vallisneri 2008; Moore Sz Hellings 



2002[|Crowder|2006[ ). This section reviews the covariance analysis for a linear least squares process, 
based on the Fisher information matrix, the method which forms the core of the error analysis in 
this paper. 

Let us suppose that the LISA combined data stream consists of discrete samples of a signal 
given by (Eq. |4|, with added noise: 



ya(Ai) + Tic 



(6) 



Here, Aj are the parameters of the source and is the noise, assumed to be stationary and 
Gaussian. The probability distribution of the ath data point is therefore 



p{sa\Xi) 



1 



a/2 



TTa„ 



(7) 



where cTq (with Greek subscript) is the standard deviation of the noise in the ath data point. 

The likelihood function for a particular data set, with parameters Aj, is the product of the 
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probabilities (Eq. [7| for each data point. It is 

L(sa| Aj) oc exp 



1 [Sa -ya(Aj)]^ 



(8) 



The set of parameters, Aj, that maximizes the hkehhood function is an unbiased estimate of the 
the set of actual model parameters, Aj. To calculate the Aj, we assume that the differences between 
the estimated values and the true values, AAj = Aj — Aj are small enough that ya{K) can be 
approximated by its first-order Taylor series expansion about A^: 



ya(Afc) « ya{Xk) + ViyoAAi 



(9) 



where Vj represents the partial derivative with respect to Aj. This first-order expansion is valid 
when the SNR is high enough, and the degree of correlation between the parameters is low enough 
that the resulting AAj are small. Using Eq. |8]we fincQthat the likelihood is maximized by 



AA,: 



where the matrix F is the Fisher information matrix^ with components 



1 



The expected parameter covariance matrix is 

(A\AA,-) 



F-' 



(10) 



(11) 



(12) 



The standard deviations in each detected parameter, cjj (with Latin subscript), are given by the 
diagonal elements of the covariance matrix: 

o\ = F^^^ (no sum over i) (13) 

It is important to remember that the Fisher error estimate is accurate only when the parameter 



uncertainties are small compared to the characteristic scales of the system being fit (Vallisneri &; 



Mock LISA Data Challenge Taskforce 20061, a condition that is not well satisfied for all of the 



binaries being modelled here. In these cases, however, the method tends to overestimate the degree 
of uncertainty in systems with a sharply-defined minimum, and the resulting error estimates tend 
to be conservative. 

Rather than write our own Fisher error estimation codes, we have made use of the publicly 



available LISA Calculator (Crowder & Cornish 2006 Crowder 20061. The LISA Calculator uses 



the same instrument noise model Sn{f) that is used as input to the SCG, and an analytic signal 



'it is also necessary to keep only terms that are first order in the inverse of the SNR; see 



Vallisneri 



(20081 
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model given by Eqs. 3 and 4. It takes as input a set of source parameter values, Xj, and outputs a 
set of standard deviations (equation 13), crj, for the detected parameter values, Aj. Each detected 
parameter, Aj, is assumed to have a Gaussian probability density with mean Aj and standard 
deviation cjj, 

l2" 



PD{Xi\K) 



1 



2TTaf 



: exp 



[A— A ,]- 
2a? 



(14) 



4. The Error Kernel 

The output of the astrophysical models of interest can be described in terms of a coalescence 
rate, T{M,ri,z), per unit redshift and time. However, the coalescence rate observed in the LISA 
detector, T'(M, fj, z), will differ from T(M, r], z) because some sources will be too weak to be detected 
and because errors in the LISA parameter determination will assign incorrect parameters to the 
source, due to the effect of the noise on the estimation process. The effect of these errors is summed 
up in the LISA error kernel, K: 

T'{M,fi,z) = f T{M,r],z)xeiM,r],z) 

J pop 

xK{M, fj, z\M, T], z) dz dri dM, (15) 

where e(M, rj, z) is the average detectability of a source with parameters {M, 77, z} in the LISA 
detector. 



4.1. Calculating the Error Kernel 

As discussed in Section 3.1, coalescence rates given by the various models are functions of the 
population parameters only. They do not depend on the sample parameters, which arise from the 
random relationship between the observer and a particular binary in the population. We therefore 
produce a Monte Carlo average or 'marginalisation' over the sample parameters, compiling the 
Fisher matrix error estimates into an 'Error Kernel' which is a function of sample parameters only. 

Since we have no a priori reason to expect inhomogeneous or anisotropic distribution, the 
values of the extrinsic sample parameters of a black hole binary are assumed to be uniformly 
distributed - angular location and orientation variables are uniformly distributed on the sky, and 
^0 is uniformly distributed over the interval [0,27r] (see Table [T|. 

The appropriate distribution to use for tc is somewhat more complicated, owing to two primary 
considerations. First, astrophysical models of the MBH population are usually expressed in terms 
of the number of coalescences per unit time and redshift. Second, the LISA detectability of a 
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binary (generally related to SNR) is a monotonically decreasing function of the binary's tc (all 
else being equal) for tc longer than the LISA observation time Tots- Because of their relatively 
stronger signals, binaries coalescing inside or soon after the LISA observation lifetime window are 
far more likely to be detected than those with long times to coalescence, so those with long times 
to coalescence represent a negligible fraction of the set of detected binaries. In calculating the error 
kernel, we found that the number of additional detected sources coalescing in the year following 
the end of LISA observation was a negligible fraction of the total, and decided to simply use a 1 
year observation time and range of tc for most of our results, multiplying by 3 to get results for 
a 3 year observation. We found that the results were not significantly different from, for instance, 
using a full 3 year observation time and 4 year range of tc. 

While the binary masses are also population parameters, the models we have found in the 
literature generally divide their mass spectra into very wide logarithmic bins or give no mass spectra 
at all. Since the mass also has a significant effect on the detectability and parameter estimation 
error of a source, however, we cannot simply generate masses completely at random or treat them 
in some other trivial fashion. We therefore partially marginalize over the mass and attempt to 
make a reasonable choice for the mass distribution within the mass bins published in the literature. 
This choice of mass distribution is somewhat problematic, since additional information on the mass 



distribution is available for some of the models studied (Volonteri et al. (20031, for instance), but not 
for others. Even if useful information on the mass distribution were available for all of the models 
studied, we prefer not to tailor intra-bin mass distributions to each particular model, because we 
want the error kernels to be model-independent. We have decided, for purposes of this paper, to 
use a simple uniform logarithmic distribution within each mass bin. While this distribution can 
produce significant differences in detection rates for the coarsely-binned models studied here, our 
opinion is that it remains a reasonable choice for a model independent intra-bin mass distribution, 
and such problems are best solved by increasing the mass resolution of the reported model results. 
For similar reasons, we also completely marginalise over mass ratio. We use the three mass ranges 



found in Sesana et al. (2007) for our mass bins. The redshift is not marginalised, and separate 
Monte Carlo runs are made at uniformly spaced values of z ranging from ... 20 (see Table [T]) 
and the statistics are collected as a function of redshift. The 'population' parameter space is thus 
considered as a collection of two-dimensional volume elements chosen from the three mass bins and 
up to 80 bins in redshift. 

Within each volume element, values for the marginalised parameters are chosen randomly and 
a covariance error analysis is performed, with the probability distribution function (PDF) for each 



sample being calculated using Eq. 14 A typical single PDF for one sample is shown in Figure 
la. PDFs from the random sampling of the source parameters are stored and averaged, producing 
a PDF for the population parameters corresponding to the source volume element chosen. An 
example of a marginalised PDF is shown in Figure lb. This conditional probability - the probability 
of LISA assigning a particular redshift to a source, given that the source parameters are within this 



mass and redshift bin - is exactly what is meant by the error kernel in Eq. 15 The marginalised 
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error kernel is thus 



N , r rt , n2 



(16) 



where (Tj^^ is the uncertainty in the ith parameter in the /xth randomly generated set of source 
parameters within the bin. N is the total number of samples generated in the bin. 

Although each individual element of the sum is a Gaussian, the process results in a non- 
Gaussian distribution since the size of the uncertainty, CTj^^, changes with each new set of sample 
parameters. This can be seen by comparing Figures [T^ andjT]^. The sum of Gaussians, all centred 
on the value of the source parameter, has a taller peak and fatter tails than a single Gaussian with 
average standard deviation [^^^i^-j?^] ' /N. 

There are two limitations in the way we have generated the error kernel that stem from our 
use of the linear least-squares LISA Calculator. First, sources at moderate redshift and large sigma 
will have tails that extend to low z, even though a true nearby source would not be confused with 
a stronger source at moderate redshift. Production of the true PDF for such a case would involve 
a more complete algorithm, avoiding the limitations of linear least-squares analysis and resulting 
in a shorter low-z tail. Second, the LISA Calculator, like many least-squares tools, drops an ill- 
determined parameter when the information matrix is singular, and gives an inappropriately low 
sigma for the remaining parameters. This did occur in a number of the cases we ran and contributed 
some anomalously strong peaks to the Monte Carlo averaging. 

Our covariance studies found that the fractional uncertainties in the redshifted mass variables 
were always much less than the fractional uncertainties in redshift (see figure|2]), and are insignificant 
compared to our coarse mass binning For the purposes of this initial study, we have ignored the 
mass errors and considered only the distribution of the detected redshifts. The error kernel is thus 



(17) 



where i corresponds to one of the mass ranges defined in Table [T] and where the sigmas are under- 
stood to be the uncertainties determined for each mass bin. Thus, our Monte Carlo study returns 
three error kernels at a given redshift, one for each mass range. Each error kernel consists of a set 
of PDFs, one PDF for each redshift bin z. These source redshifts are uniformly spaced between 
and 20 (see figure [TJ^). Each of these PDFs represent Monte Carlo averages of LISA calculator 
Gaussian PDFs, varied over the sample parameters and the mass parameter ranges listed in Table 
[T] Each error kernel contains about 2 million LISA calculator runs. One may get a feel for the 
shape of the entire kernel by looking at the 70% confidence intervals shown in Figure [l]l. 



''Rest-frame variables with dimensions of mass will have an error that is 100% correlated with the error in the 
redshift, a detail we will investigate in future work. 
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source z 

Fig. 1. — Creating the error kernel for the 'medium mass' bin and marginalised parameters described 
in Table [l] total masses between lOjOOOM© and IOO.OOOMq. a) The Gaussian PDF imphed by 
simple RMS error, b) Adding the probability densities resulting from many different marginalised 
parameters results in a highly non-normal distribution. Shown here (in solid black) is the dis- 
tribution of possible detections given several hundred sources at a redshift of 10 in a mass range 
10^ : IO^Mq. Overlaid (in dashes) is the distribution obtained by simply adding the errors in 
quadrature, c) The Zs = 10 PDF inserted into its place in the error kernel. Source redshifts are 
sampled at even redshift intervals of 0.25. d) 70% confidence intervals for LISA determination of 
redshift gives an overview of the resulting error kernel. 
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Fig. 2. — Average fractional error for binaries witli total mass between 10, OOOMq and 1, 000, OOOM© 
and reduced mass between 0.001 and 0.25, using + for redshift and x for reduced mass. Chirp mass 
errors are not shown, but are an order of magnitude below the reduced mass errors. Error in 
redshift dominates the reduced mass error by over an order of magnitude. 
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4.2. Applying the Error Kernel 



Once the error kernel in Eq. 17 has been calculated, it may be applied to any population model 
that gives the source coalescence rate in the corresponding mass {i) and redshift bin Ti{z), producing 
a prediction for the detected coalescence rate, T'-(z). This is accomplished by a straightforward 
convolution, which we write in continuum form as 

T[{z) = / T,{z)xKi{z,z)dz. (18) 



As we noted in Eq. [T5| some of the binaries sampled will have signal-to-noise ratios too small 
to be detectable, regardless of the error in the parameters. Previous analysis by [Sesana et al. 



(20071 has used the fiducial SNR limit of 5. We chose to use a cutoff at SNR = 8, but the number 
of additional sources dropped due to our more conservative cutoff was negligible. The PDFs of 
binaries which do not make the SNR cutoff should not be included in the error kernels, but the 
proportion of rejected binaries in each source parameter bin must still be taken into account when 
the error kernels are applied to the models. Therefore, the error kernels (which would otherwise 
normalise to 1) are weighted with the fraction of binaries in each bin, ei(z), having SNR > 8. 
Taking this detectability into account, the convolved coalescence rates may be written 

T[{z)= Ti{z)xei{z)x Ki{z,z)dz. (19) 



5. Discriminating Between Population Models 



As an illustration of LISA's ability to discriminate between black hole population models, 
we consider four formation models. The models we have chosen for the demonstration - those by 
Volonteri et al. (20031, Koushiappas et al. ( 2004[ ), and two bylBegelman et al. (20061, one with 'high' 
feedback and one with 'low' feedback, (hereafter VHM, KBD, BVRhf, and BVRlf respectively) - 
are variations on the extended Press-Schecter (EPS) formalism by Lacey &: Cole ( 1993| ) which 
assigns a mass-dependent probability to halo mergers. Key variations between these models are 
their assumptions for accretion, their binary hardening scenarios, their choices for mass and redshift 
of seed formation, and the details of the way they handle MBH binary interactions near the merger. 



5.1. Convolving The Models w\t\\ the LISA Error Kernel 

The effect of the error kernel on population model testing can be seen in Figure |3] where we 
have taken the VHM model with its three mass bins lumped together, spanning the range from 
300 Mq to 10*^ M0. The VHM model, with its unique seeding scenario, predicts a large number 
of low mass, high redshift binaries. The distribution, shown as the large-amplitude solid curve in 
Figure |3] peaks and then rolls off sharply at z ~ 17. The number of sources in this model that are 
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expected to be visible with LISA, using a SNR> 8 cut-off relative to the averaged sensitivity curve 



from the SCG, is given by the low-amplitude dotted curve (see Sesana et al. 2007). The dashed 
curve represents our results, produced by integrating the model with the kernel, as in Eq. 18 It 
gives predictions for the distribution of best-fitting parameters detected by LISA, using a cutoff 
at SNR = 8. The obvious point to be made is that the redshift uncertainties smear out a model's 
features, so that, while the VHM model of Figure |3] has a very distinctive shape in its theoretical 
incarnation, the distribution that would be observed by LISA will be far less so. 

The presence of detected sources at very low {z < 2) redshifts, where the VHM population 
model says that there should be very few, is a result of the excessively large low-z tails of the error 
kernel discussed in Section |4| The less sharply peaked shape of the error kernel results as compared 
with the Sesana results, however, is a consquence of applying detection errors to the redshifts, and 
will persist even when more robust error-estimation techniques are employed. 



Even with our SNR cutoff, we predict more visible sources (179 vs. 96) than did Sesana et al 



( 2007 1 . This is due in part to our error kernel, with its use of the logarithmic mass distributions 



within the large mass bins found in the literature, having more massive binaries at high redshift 



than is actually the case for the VHM models. We also found a puzzling discrepancy in Sesana 



et al. (20071 between their stated event counts ( 250 for the VHM model with a 3 year range of 
coalescence times) and the event counts found by integrating the curves in their figure 1 (over 400 
for the same case). Since our model event rates were obtained by extracting the curves from that 
figure, this discrepancy could also contribute to the differing number of visible sources. 



5.2. Discriminating Between Models 

For the four models we have chosen to consider as illustrative examples, the results of the 
error kernel convolutions are shown in Figure |4j The graph in the upper left is for all masses and 
the other three graphs represent the three mass bins we used. We use a modified version of the 
Kolmogorov-Smirnov (K-S) test as a measure of separability of the models. Our test differs from 
the K-S test in that it is sensitive to differences in the model event rates as well as to the cumulative 
distribution functions (CDFs) of samples drawn from the models. For each pair of models shown in 
Table |2] we have simulated Monte Carlo draws of the number of sources in each redshift bin, using a 
Poisson distribution with probability given by each of the models. One thousand draws were taken 
for each model and our test statistic was calculated, finding the greatest deviation between the 
cumulative histograms of the two models. The probability that the two draws were from the same 
model was then found by Monte Carlo sampling from the null hypothesis that the two samples have 
the same CDF and have event counts which are Poisson distributed with identical rate parameters. 



Several comparisons were done between the four models chosen from Sesana et al. (2007 1 . 
For each comparison, we assumed one year of LISA observations and randomly drew coalescence 
parameters using the probability distributions for the two models being compared. In each set of 
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Raw VHM Event Counts 

Detection counts: SNR cut and error kernel 




redshift 



Fig. 3. — The effect of the LISA error kernel on the MBH binary population predicted by Volonteri 



et al. (20031 for a 3 year observation. Solid line: modelled source distribution for masses below 



IO^Mq. Dotted line: modelled source distribution with SNR > 8 cut, as applied by Sesana et al. 



(2007 1 . Dashed line: convolution of modelled source distribution with LISA error kernel. The error 
kernel distribution incorporates the large errors inherent to binary redshift determination. 
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Fig. 4.- 
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draws, our modified K-S statistic, E, was determined and Q, the probability that two such draws 
would be produced by the same model, was calculated. The results for each random draw were then 
averaged over one thousand such realisations, giving the values displayed in Table |2] In the two 
data columns, we show the results from the raw models themselves, with no parameter uncertainties 
taken into account, and the results from the models after they have been convolved with the LISA 
error kernel. As can be seen in the table, the probability that any of the simulated data sets for 
one model might have been produced by one of the other models is small. The models examined 
here appear to be easily distinguishable from each other, with the exception of the comparison of 
the BVRhf model with the BVRlf model with its average Q value of 0.055 (corresponding to a 
rather shaky 94.5% confidence). Even in that case, when we look at the median of Q rather than 
its mean, we find it to be 0.012 (98.8% confidence), implying that a BVRlf realization can usually 
be distinguished from a BVRhf realization. 



6. Conclusion 

In this paper, we have constructed an 'error kernel' for the LISA detector, using a public- 
domain error computation module in a Monte-Carlo data pipeline, marginalising over the 'sample' 
parameters. This result represents a first implementation of a model-level error function for a 
gravitational wave observatory, a concept which is central to all types of astronomy. The error 
kernel approach introduced here is designed to replace simple SNR cuts as the interface between 
population modellers and gravitational wave signal specialists. 

The error kernel can be computed as a function of any of the population parameters. We 
have chosen to use redshift alone, but total mass or mass ratio could easily have been added. 
Once calculated, contours of the error kernel can be used to visualise the relative impact of the 
detector noise, analysis technique, and choice of parameter binning on the ability of the detector 
to determine the parameters of the sources in the model. 

We reiterate that we are limited by working with the models as published, with their emphasis 
on redshift as the parameter of interest, rather than having access to more detailed model results 
giving populations as functions of masses and redshifts. In particular, while the logarithmic intra- 
bin mass distribution is as reasonable a choice as any other, the mass bins published in the literature 
are so large that no model-independent distribution can produce an error kernel that is free from 
significant bias compared to what would be obtained using the full model results. Furthermore, 
with LISA'S exquisite resolution in the redshifted mass variables, the mass distribution is likely to 
provide useful astrophysical information in its own right. It is our opinion that a mass resolution 
of at least 1 bin per decade over the range 10^ . . . 10^ Mq is necessary to meaningfully specify the 
mass dependence of these populations. We are currently working toward improved comparisons 
using more detailed models. 

The demonstration in this paper for four model comparisons looks at the distinguishability of 
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Binary Parameter 


Marginalisation Range(s) 






: vr 




L 


: 7r/2 







: 27r 






-7r/2 : 7r/2 




$0 


: vr 










■n 


0.0025 : 0.25 




z 


Not Marginalised 






250 : lO^Mo 


Low Mass Case 


Mtot 


10^ : IO^Mq 


Medium Mass Case 




10^ : lO^Mo 


High Mass Case 



Table 1: This table of parameters lists the completely marginalised case parameters (grey boxes) 
along with their ranges and describes the treatment of the three population parameters. Reduced 
mass ratio 77 is averaged, z is held constant and total mass {Mtot) is divided into three units. 



Models Compared 


Before LISA kernel 


After LISA Kernel 






{Q) 


{E) 


{Q) 


VHM - KBD 


175.6 


< 10-4 


90.5 


< 10-4 


VHM - BVRlf 


119.7 


< 10-4 


49.7 


< 10-4 


VHM - BVRhf 


132.9 


< 10-4 


58.9 


< 1.0 X 10-4 


BVRhf - BVRlf 


14.29 


0.021 


9.82 


0.055 



Table 2: Comparisons of models before and after convolution with LISA kernel, for binaries coa- 
lescing within the observation window, assuming one year of observation time. The 'Before LISA' 
comparisons effectively assume that all of the sources are detectable and have zero redshift error, 
while the 'After LISA' comparisons incorporates the effects of both parameter uncertainty and de- 
tectability. E is the maximum deviation between the cumulative histograms of random draws from 
the two models. Q is the corresponding probability that random fluctuations could be responsible 
for the deviation. 
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the models based on what might eventually be a catalogue of LISA binaries. Our approach has been 
to do a forward modelling from the relevant parameters of the population model to the detected 
parameters of sources seen in the LISA data. The fact that those predictions are statistically 
different suggests that LISA data will have appreciable model-discriminating power. 

When the LISA data are finally available, the analysis will include backward, Bayesian, mod- 
elling which calculates and interprets the contribution of each detected source to a likelihood func- 
tion for the models being tested. This backward population analysis framework will be an essential 
component of future LISA data analysis, and is a natural direction of our future work in the long 
term. 
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